#!/usr/bin/env python
# =========== Enter 4 corner's coordinates(Basic Nodes) ======================
# in clockwise or counter colockwise order
# Basic node 'a'
xa=3
ya=-0.1
za=0.1

# Basic node 'b'
xb=0
yb=-0.1
zb=0.1

# Basic node 'c'
xc=0
yc=0.1
zc=0.1

# Basic node 'd'
xd=3
yd=0.1
zd=0.1

# =========== Enter dividing numbers ==============================
# number of elements between 1st and 2nd nodes
n1=90
# number of elements between 2nd and 3rd nodes
n2=6

# =========== Enter starting number of nodes and elements =============
sn=7060
se=1621

# =========== Enter file path and File name =======================================

Path='D:\\PyProgram'
Filename='Mesh4.tcl'

# =========== Type of Shell Element ============================
TypeofEle='ShellNL'
SectionNum='4'

# ============== Input complete ===================================

# =========== Inserting nodes into the 4 boundaries =======================

Nab = [[] for i in range(2*n1+1)]
Nbc = [[] for i in range(2*n2+1)]
Ncd = [[] for i in range(2*n1+1)]
Nda = [[] for i in range(2*n2+1)]

for i in range(0,2*n1+1):
    x1=xa+float(xb-xa)/(2*n1)*i
    y1=ya+float(yb-ya)/(2*n1)*i
    z1=za+float(zb-za)/(2*n1)*i
    Nab[i].append(x1)
    Nab[i].append(y1)
    Nab[i].append(z1)
    
for i in range(0,2*n2+1):
    x1=xb+float(xc-xb)/(2*n2)*i
    y1=yb+float(yc-yb)/(2*n2)*i
    z1=zb+float(zc-zb)/(2*n2)*i
    Nbc[i].append(x1)
    Nbc[i].append(y1)
    Nbc[i].append(z1)

for i in range(0,2*n1+1):
    x1=xc+float(xd-xc)/(2*n1)*i
    y1=yc+float(yd-yc)/(2*n1)*i
    z1=zc+float(zd-zc)/(2*n1)*i
    Ncd[i].append(x1)
    Ncd[i].append(y1)
    Ncd[i].append(z1)
    
for i in range(0,2*n2+1):
    x1=xd+float(xa-xd)/(2*n2)*i
    y1=yd+float(ya-yd)/(2*n2)*i
    z1=zd+float(za-zd)/(2*n2)*i
    Nda[i].append(x1)
    Nda[i].append(y1)
    Nda[i].append(z1)
    
#============ Inserting nodes into area =============

Narea = [[] for i in range((2*n1+1)*(2*n2+1))]
for i in range(0,2*n1+1):
    Narea[i].append(Nab[i][0])
    Narea[i].append(Nab[i][1])
    Narea[i].append(Nab[i][2])
    
for j in range(1,2*n2):
    for i in range(0,2*n1+1):
        cg=Nab[i][0]
        cg1=Ncd[2*n1-i][0]
        cg3=Nbc[j][0]
        cg5=Nda[2*n2-j][0]
        cg7=Nab[i][1]
        cg9=Ncd[2*n1-i][1]
        cg11=Nbc[j][1]
        cg13=Nda[2*n2-j][1]
        t1 = cg * cg3
        t4 = cg * cg5
        t7 = cg1 * cg3
        t10 = cg1 * cg5
        if (cg * cg11 - cg * cg13 - cg1 * cg11 + cg1 * cg13 - cg3 * cg7 + cg3 * cg9 + cg5 * cg7 - cg5 * cg9)!=0:
            x1 = (t1 * cg9 - t1 * cg13 + t10 * cg7 - t10 * cg11 - t4 * cg9 + t4 * cg11 - t7 * cg7 + t7 * cg13) / (cg * cg11 - cg * cg13 - cg1 * cg11 + cg1 * cg13 - cg3 * cg7 + cg3 * cg9 + cg5 * cg7 - cg5 * cg9)
        
            tt1 = cg * cg9
            tt4 = cg1 * cg7
            tt7 = cg3 * cg7
            tt9 = cg3 * cg9
            tt11 = cg5 * cg7
            tt13 = cg5 * cg9
            y1 = (tt1 * cg11 - tt1 * cg13 + tt11 * cg11 - tt13 * cg11 - tt4 * cg11 + tt4 * cg13 - tt7 * cg13 + tt9 * cg13) / (cg * cg11 - cg * cg13 - cg1 * cg11 + cg1 * cg13 + tt11 - tt13 - tt7 + tt9)
            if(cg1-cg!=0):
                z1=Nab[i][2]+(Ncd[i][2]-Nab[i][2])*float(x1-cg)/(cg1-cg)
            elif(cg5-cg3!=0):
                z1=Nab[i][2]+(Ncd[i][2]-Nab[i][2])*float(y1-cg3)/(cg5-cg3)
        
        else:
            cg=Nab[i][1]
            cg1=Ncd[2*n1-i][1]
            cg3=Nbc[j][1]
            cg5=Nda[2*n2-j][1]
            cg7=Nab[i][2]
            cg9=Ncd[2*n1-i][2]
            cg11=Nbc[j][2]
            cg13=Nda[2*n2-j][2]
            t1 = cg * cg3
            t4 = cg * cg5
            t7 = cg1 * cg3
            t10 = cg1 * cg5
            if (cg * cg11 - cg * cg13 - cg1 * cg11 + cg1 * cg13 - cg3 * cg7 + cg3 * cg9 + cg5 * cg7 - cg5 * cg9)!=0:
                y1 = (t1 * cg9 - t1 * cg13 + t10 * cg7 - t10 * cg11 - t4 * cg9 + t4 * cg11 - t7 * cg7 + t7 * cg13) / (cg * cg11 - cg * cg13 - cg1 * cg11 + cg1 * cg13 - cg3 * cg7 + cg3 * cg9 + cg5 * cg7 - cg5 * cg9)
        
                tt1 = cg * cg9
                tt4 = cg1 * cg7
                tt7 = cg3 * cg7
                tt9 = cg3 * cg9
                tt11 = cg5 * cg7
                tt13 = cg5 * cg9
                z1 = (tt1 * cg11 - tt1 * cg13 + tt11 * cg11 - tt13 * cg11 - tt4 * cg11 + tt4 * cg13 - tt7 * cg13 + tt9 * cg13) / (cg * cg11 - cg * cg13 - cg1 * cg11 + cg1 * cg13 + tt11 - tt13 - tt7 + tt9)
                if(cg1-cg!=0):
                    x1=Nab[i][0]+(Ncd[i][0]-Nab[i][0])*float(x1-cg)/(cg1-cg)
                elif(cg5-cg3!=0):
                    x1=Nab[i][0]+(Ncd[i][0]-Nab[i][0])*float(y1-cg3)/(cg5-cg3)
            
            else:
                cg=Nab[i][2]
                cg1=Ncd[2*n1-i][2]
                cg3=Nbc[j][2]
                cg5=Nda[2*n2-j][2]
                cg7=Nab[i][0]
                cg9=Ncd[2*n1-i][0]
                cg11=Nbc[j][0]
                cg13=Nda[2*n2-j][0]
                t1 = cg * cg3
                t4 = cg * cg5
                t7 = cg1 * cg3
                t10 = cg1 * cg5
                z1 = (t1 * cg9 - t1 * cg13 + t10 * cg7 - t10 * cg11 - t4 * cg9 + t4 * cg11 - t7 * cg7 + t7 * cg13) / (cg * cg11 - cg * cg13 - cg1 * cg11 + cg1 * cg13 - cg3 * cg7 + cg3 * cg9 + cg5 * cg7 - cg5 * cg9)
                tt1 = cg * cg9
                tt4 = cg1 * cg7
                tt7 = cg3 * cg7
                tt9 = cg3 * cg9
                tt11 = cg5 * cg7
                tt13 = cg5 * cg9
                x1 = (tt1 * cg11 - tt1 * cg13 + tt11 * cg11 - tt13 * cg11 - tt4 * cg11 + tt4 * cg13 - tt7 * cg13 + tt9 * cg13) / (cg * cg11 - cg * cg13 - cg1 * cg11 + cg1 * cg13 + tt11 - tt13 - tt7 + tt9)
                if(cg1-cg!=0):
                    y1=Nab[i][1]+(Ncd[i][1]-Nab[i][1])*float(x1-cg)/(cg1-cg)
                elif(cg5-cg3!=0):
                    y1=Nab[i][1]+(Ncd[i][1]-Nab[i][1])*float(y1-cg3)/(cg5-cg3)
        
        Narea[(2*n1+1)*j+i].append(x1)
        Narea[(2*n1+1)*j+i].append(y1)
        Narea[(2*n1+1)*j+i].append(z1)

j=2*n1

for i in range(-2*n1-1,0):
    Narea[i].append(Ncd[j][0])
    Narea[i].append(Ncd[j][1])
    Narea[i].append(Ncd[j][2])
    j=j-1
    
print(Narea)

#======================== Create Elements ==============================

Elements=[]
for i in range(0,(2*n1+1)*(2*n2-1),2):
    if i%(4*n1+2)<=2*n1-2:
        ele=[int(i),int(i+2),int(i+4+4*n1),int(i+2+4*n1),int(i+1),int(i+2*n1+3),int(i+3+4*n1),int(i+2*n1+1),int(i+2*n1+2)]
        Elements.append(ele)
        
print(Elements)

# ====================== Writing nodes to tcl file ======================

fn=Path+'/'+Filename
fn=open(fn,'w')
i=sn
for obj in Narea:
    fn.write('node %d %f %f %f'%(i,obj[0],obj[1],obj[2]))
    fn.write('\n')
    i=i+1
   
i=se 
for obj in Elements:
    fn.write('element %s %d %d %d %d %d %d %d %d %d %d %s'%(TypeofEle,i,obj[0]+sn,obj[1]+sn,obj[2]+sn,obj[3]+sn,obj[4]+sn,obj[5]+sn,obj[6]+sn,obj[7]+sn,obj[8]+sn,SectionNum))
    fn.write('\n')
    i=i+1
    
fn.close()
    